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Abstract. The hybrid kinetic model supports comprehensive simulation 
of the interaction between different spatial and energetic elements of the Europa 
moon-magnetosphere system with respect a to variable upstream magnetic field and 
flux or density distributions of plasma and energetic ions, electrons, and neutral atoms. 
This capability is critical for improving the interpretation of the existing Europa flyby 
measurements from the Galileo Orbiter mission, and for planning flyby and orbital 
measurements (including the surface and atmospheric compositions) for future missions. 
The simulations are based on recent models of the atmosphere of Europa (Cassidy et 
al., 2007; Shematovich et al., 2005). In contrast to previous approaches with MHD 
simulations, the hybrid model allows us to fully take into account the finite gyroradius 
effect and electron pressure, and to correctly estimate the ion velocity distribution and 
the fiuxes along the magnetic field (assuming an initial Maxwellian velocity distribution 
for upstream background ions). Photoionization, electron-impact ionization, charge 
exchange and collisions between the ions and neutrals are also included in our model. 
We consider the models with O"*""*" and S~^~^ background plasma, and various betas for 
background ions and electrons, and pickup electrons. The majority of O2 atmosphere is 
thermal with an extended non-thermal population (Cassidy et al., 2007). In this paper 
we discuss two tasks: (1) the plasma wake structure dependence on the parameters 
of the upstream plasma and Europa's atmosphere (model I, cases (a) and (b) with a 
homogeneous Jovian magnetosphere field, an inductive magnetic dipole and high oceanic 
shell conductivity); and (2) estimation of the possible effect of an induced magnetic 
field arising from oceanic shell conductivity. This effect was estimated based on the 
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difference between tlie observed and modeled magnetic fields (model II, case (c) with 
an inhomogeneous Jovian magnetosphere field, an inductive magnetic dipole and low 
oceanic shell conductivity). 

Keywords: Europa, Jovian magnetosphere. Plasma, Magnetic fields. Ion 
composition 
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1. Introduction 

The interaction of the Jovian plasma torus with Europa and other moons is a 

fundamental problem in magnetospheric physics (see e.g., Goertz, 1980; Southwood 
et al., 1980; Southwood et al., 1984; Wolf-Gladrow et al., 1987; Ip, 1990; Schreier et 
al., 1993; Lellouch, 1996). The plasma environment near Europa was studied by flyby 
observations during the Galileo prime mission and the extended Galileo Europa mission 
(Kivelson et al., 1997; Khurana et al., 1998; Kivelson et al., 1999, Paterson et al., 1999). 

Europa, one of the icy moons of Jupiter, was encountered by the Galileo satellite 
three times during its primary mission, seven times during its Galileo Europa Mission 
(GEM), and once during Galileo Millennium Mission (GMM). Europa is located at a 
radial distance of 9.4 Rj (Jovian radii, 71,492 km) from Jupiter, and has a radius of 
1560km (1 Re). 

The interaction of Europa with the magnetized plasma of the Jovian plasma sheet 
gives rise to a so-called Alfven wing, which has been extensively studied in the case of 
lo (e.g., Neubauer, 1980; Southwood et al., 1980; Herbert, 1985; Lipatov and Combi, 
2006). Neubauer (1998; 1999) has shown theoretically how an Alfven wing is modified 
by an induced magnetic field, such as that found at Europa (Kivelson et al., 2000). 
Observations by Kivelson et al. (1992) show the generation of ultra-low frequency 
electromagnetic waves in Europa's wake. These waves have frequencies near and below 
the gyrofrequencies of the ion species in the plasma torus (e.g., ionized sulfur, oxygen, 
and protons). Ion cyclotron waves grow when ion distribution functions are sufficiently 
anisotropic, as occurs when ion pickup creates a ring distribution of ions (in velocity 
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space). The analysis of these waves has been done by Huddleston et al. (1997) (lo), 
Volwerk et al. (2001) and Kivelson, Khurana and Volwerk (2009) (Europa). They found 
intensive wave power at low frequencies (near and below the cyclotron frequencies of 
heavy ions) in Europa's wake during the Ell and E15 fiybys. However, our current 
3D hybrid modeling cannot yet produce these waves due to insufficient spatial grid 
resolution. 

The most general and accurate theoretical approach to this problem would require 
the solution of a nonlinear coupled set of integro-MHD/kinetic-Boltzmann equations 
which describe the dynamics of Jupiter's corotating magnetospheric plasma, pickup 
ions, and ionosphere, together with the neutrals from Europa's atmosphere. To first 
order, the plasma and neutral atoms and molecules are coupled by charge exchange and 
ionization. The characteristic scale of the ionized components is usually determined 
by the typical ion gyroradius, which for Europa is much less than characteristic global 
magnetospheric scales of interest, but which may be comparable to the thickness of 
the plasma structures near Europa. Kinetic approaches, such as Direct Simulation 
Monte Carlo, have been apphed to the understanding of global aspects of the neutral 
atmosphere (Marconi et al., 1996; Austin and Goldstein, 2000). Plasma kinetic modeling 
is, however, much more complicated, and even at the current stage of computational 
technology require some approximations and compromises to make some initial progress. 
Several approaches have been formulated for including the neutral component and 
pickup ions self-consistently in models that describe the interaction of the plasma torus 
with Europa. 
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There have been recent efforts to improve and extend the pre-Gahleo models for 
Europa, lo and Ganymede, both in terms of the MHD (Kabin et al., 1999; Combi et al., 
1998; Linker et al., 1998; Kabin et al., 2001; Jia et al., 2008), the electrodynamic (Saur 
et al., 1998; Saur et al., 1999; Schilling et al., 2008), and hybrid kinetic (Lipatov and 
Combi, 2006; Lipatov et al., 2010) approaches. These approaches are distinguished by 
the physical assumptions that they include. MHD and hybrid kinetic models cannot, 
at least yet, include the charge separation effects which are likely to be important very 
close to the moon where the neutral densities are large and the electric potential can 
introduce non-symmetric flow around the body. MHD models for lo either include 
constant artificial conductivity (Linker et., 1998) or assume perfect conductivity (Combi 
et al., 1998). Comparisons of the sets of published results do not indicate that this 
choice has any important consequences. The MHD model of Europa developed by Kabin 
et al. (1999) includes an cxospheric mass loading, ion-neutral charge exchange, and 
recombination. Further development of this model by Liu et al. (2000) already includes 
a possible intrinsic dipole magnetic field of Europa. Schilling, Neubauer and Saur (2007; 
2008) found that a conductivity of Europa's ocean of 500 mS/m or largecombined with 
an ocean thickness of 100 km or smaller is most suitable for explaining the magnetic 
fiyby data. They also found that the influence of the fields induced by the time variable 
plasma interaction is small compared to the induction caused by the time-varying 
background field. 

Hybrid kinetic models can include the finite ion gyroradius effects, non-Maxwellian 
velocity distribution for ions, and correct flux of pickup ions along the magnetic fleld. 
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Hybrid modeling of lo has demonstrated several features. The kinetic behavior of ion 
dynamics reproduces the inverse structure of the magnetic field (due to drift current) 
which cannot be explained by standard MHD or electro dynamic modeling which do not 
account for anisotropic ion pressure. The diamagnetic effect of non-isotropic gyrating 
pickup ions broadens the B-field perturbation and produces increased temperatures in 
the flanks of the wake, as observed by the Galileo spacecraft, but had not been explained 
by previous models. The temperatures of the electrons which are created and cooled 
by collisions with neutrals in the exosphere and inside the ionosphere may strongly 
affect the pickup ion dynamics along the magnetic fleld and consequently the pickup 
distribution across the wake. The physical chemistry in Jo's corona was considered in 
the paper by Dols et al. (2008). They couple a model of the plasma flow around lo 
with a multi-species chemistry model and compare the model results to the Galileo 
observation in lo's wake. 

Gahleo flyby measurements E4, E6 (plasma only). Ell, E12, E14, E15, E19, and 
E26 demonstrate several features in the plasma environment: Alfven wing formation and 
an induced magnetosphere, possible existence of the dipole-type induced magnetic fleld, 
and variation of the magnetic fleld in the plasma wake due to diamagnetic currents. 
The measurements also demonstrate mass loading of the plasma torus plasma by pickup 
ions and the interaction of the ions with the surface of Europa. For an interpretation of 
these data we need to use a kinetic model because of effects of the flnite ion gyroradius. 

Hybrid models have been shown to be very useful in studying the complex plasma 
wave processes of space, astrophysical, and laboratory plasmas. These models provide a 
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kinetic description of plasmas in local regions, together with the possibility of performing 
global modeling of the whole plasma system. Revolutionary advances in computational 
speed and memory are making hybrid modeling of various space plasma problems a 
much more effective general tool. 

In this paper, we apply a time-dependent Boltzmann equation (a "particle in cell" 
approach) together with a hybrid kinetic plasma (ion kinetic) model in three spatial 
dimensions (see, e.g. Lipatov and Combi, 2006; Lipatov et al., 2010), using a prescribed 
but adjustable neutral atmosphere model for Europa. A Boltzmann simulation is 
applied to model charge exchange between incoming and pickup ions and the immobile 
atmospheric neutrals. In this paper we discuss the results of the hybrid kinetic modeling 
of Europa's environment - namely, the global plasma structures (formation of the 
magnetic barrier, Alfven wing, pickup ion tail etc.). The results of these kinetic 
modeling are compared with the Galileo E4 flyby observational data. Currently, we are 
working on the hybrid model of the E12 flyby. The remarkable aspect of this flyby is a 
strong variation in the upstream plasma density profile approximately from 400 cm~^ to 
80cm~^. The results of this modeling will be discussed in future publications. 

The paper is organized as follows: in Section 2 we present the computational model 
and a formulation of the problem. In Section 3 we present the results of the modeling 
of the plasma environment near Europa and the comparison with observational data. 
Finally, in Section 4 we summarize our results and discuss the future development of 
our computational model. 
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2. Formulation of the Problem and Mathematical Model 

To study the interaction of the plasma torus with the ionized and neutral 

components of Europa's environment, we use a quasineutral hybrid model for ions 
and electrons. The model includes ionization (which in the Europa environment is 
dominated by electron impact ionization, not photoionization) and charge exchange. 
The atmosphere is considered to be an immobile component in this paper. 

In our hybrid modeling, the dynamics of upstream ions and implanted ions are 
described in a kinetic approach, while the dynamics of the electrons are described in 
a hydrodynamical approximation. The details of this plasma-neutral approach were 
developed early for the study of the lo- Jovian plasma interaction (Lipatov and Combi, 
2006). 

The single ion particle motion is described by the equations (see, e.g. Eqs. (1) and 
(14) from Mankofsky, Sudan and Denavit (1987)): 

dt dt Mi \ c J Mi Micrii 

Here we assume that the charge state is — 1. Uj, and J denote the charge- averaged 
velocity of all (incoming and pickup) ions and the total current, Eq. (5). The subscript 
s denotes the ion population (s = 1, 2 for incoming ions and s = 3, 4 for pickup ions) 
and the index / is the particle index. Uif. and are collision frequencies between ions 
and electrons, and ions and neutrals that may include Coulomb collisions and collisions 
due to particle-wave interaction. 

For a plasma, the thermal velocity, v'^ {a = i,e), is assumed greater than the drift 
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velocity, so we take 



(2) 



where the cross section cr"'" is typically about 5 x 10~^^cm^ (see, e.g., Eq. (17) from 
Mankofsky, Sudan and Denavit (1987)). 

For massless electrons the equation of motion of the electron fluid takes the form of 
the standard generalized Ohm's law (e.g. Braginskii, 1965): 



E = ^(JeXB)-— Vpe 



nip 



eUeC 



Yl ^e,s[(Uj - Us) ] + i^a,eoUe 

„ TIG 



(3) 



where Pe — nme{v'^) /?> — UekBT^, and v'^ are the scalar electron pressure and the thermal 
velocity of electrons, and the electron current is estimated from Eq. ( 5). 
The induction equation (Faraday's law) has a form 



IdB 
c~dt 



+ V X E = 0. 



(4) 



The total current is given by 



J = Je + Ji; Ji = X] e^-sUs = eUiUi, 



(5) 



where Us is the bulk velocity of ions of the type s. 



Since we suppose that electron heating due to collisions with ions is very small, the 
electron fluid is considered adiabatic. For simplicity we assume that the total electron 
pressure may be represented as a sum of partial pressures of all electron populations: 



i^e-n-l^p + /3e,pinf pi) 



(6) 
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where and ^e,Pi denote electron upwind, and pickup betas. Note that ^e,fc = 
Pe^kl {B'^ where /c is a population of electrons. We also assume here that 

^e,up — ^i,up5 ''^e.PI = "^i.PI- 

The neutral atmosphere of Europa serves as a source of new ions, mainly by 
electron impact ionization from corotating (or nearly corotating) plasma and also by 
photoionization. The neutral atmospheric molecules also serve as coUisional targets for 
charge exchange by corotating ions. The impacting ions consist of both upstream torus 
ions and newly implanted ions which are picked up by the motional electric field. 

In the current model we assume that the background plasma contains only ions 
with molecular mass/charge of 8 and 16 corresponding to O"*""*" and S^^ , respectively. 

We assume that Europa has a radius Re — 1560km. We have also adopted a 
two-species description for the neutral O2 exosphere of exponential form (Shematovich 
et al., 2005) 

'^neutraljk ~ ^atmos,k GXp 

'"exobase.k 

)/h 

atmos.k ], (7) 

where riatmos.k denotes the maximum value of the neutral density extrapolated to 
the exobase (natmos.i = 3 x lO^cm'^; natmos,2 = 8.5 x lO^cm'^; rexobase.i ~ 1700km; 
rexobase,2 ~ 1560 km), and index k denotes either non-thermal {k — 1) or thermal {k — 2) 
species. Here the scale heights /iatmos,i — 200 km and /iatmos,2 — 30 km. 

The production rate of new ions from the exosphere near Europa corresponds to 

G^exOjk OC I^i,A;^atmos,k 6Xp[ (r ^exobase,k) /^atmos,k] ) (8) 

where natmos,k denotes the value of the neutral component density at r = rexobase,k and 
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i/i^k is the effective ionization rate per atom or molecule of species k. Vi^k includes the 
photoionozation Vp^, and the electron impact ionization by the magnetospheric electrons 
^e,im- We assume that our model of the atmosphere mainly consists of O2, and we use 
the effective photoionization rate 1.7 x 10~^s~^ (Johnson et al., 2009). We also adopt 
the effective electron impact ionization rates of 2.4 x 10~®cm^/s (for 20 eV electrons) 
and 1.1 X 10~^cm^/s (for 250 eV electrons) (see e.g. Johnson et al., 2009). Since the hot 
electrons represent only 5% of the total electron density (see Voyager 1 plasma science 
(PLS) measurements analyzed by Sittler and Strobel (1987) and Bagenal (1994)) we use 
the same composition for computing the impact ionization rate. We assume that the 
Sun is located in the direction opposite the x axis. 

The interaction of ions with neutral particles by charge exchange (see Eqs. (12) - 
(15) from Lipatov and Combi, 2006) currently includes for the following reactions: 

,5++ + O2 ^ -5+ + 0+ 

+ O2 ^ O2 + Ot (9) 

The effective cross section for charge exchange {ac,ex — 2.6 x 10~^^ m^) was the same as 
that used in the hybrid modeling of lo's plasma environment (see Lipatov and Combi, 
2006; and McGrath and Johnson, 1989). A more complete hst of reactions will be 
considered in future modeling. Of course, this also requires the addition of Monte Carlo 
computations. However, this approach is beyond the scope of this paper. 
Our code solves equations (1) - (9). 
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We discuss two models of the interaction between the Jovian magnetosphere and 
Europa. In Sect. 3.1 we discuss the interaction model for the cases with different ion 
and electron betas, different pickup ion production rates near the surface of Europa, 
and a homogeneous global Jovian magnetic field (model I, cases (a) and (b)). In in 
Sect. 3.2 we consider model II, case (c) with a realistic global Jovian magnetic field 
and the internal dipole magnetic field placed in the center of Europa. To study the 
interaction of the plasma torus with the ionosphere of Europa, the following Jovian 
plasma torus and ionosphere parameters were adopted in accordance with the Galileo 
Europa E4 flyby observational data (Paterson, Prank and Ackerson, 1999; Khurana et 
al., 1998; Kivelson et al., 1997; Kivelson et al., 1998): magnetic field. Bo = 469 nT 
and B = (77.6, —140.7, — 441.3) nT; torus plasma speed relative to Europa (Paterson, 
Frank and Ackerson, 1999), Uq — 105km/s; upstream ion densities, po++ — lOcm"^; 
Ps++ = 10 cm~^ and ion temperature, Tj = (25 — 100) eV (Paterson, Frank and 
Ackerson, 1999); electron temperature for suprathermal population, T^, — 20 eV (Sittler 
and Strobel, 1987); ratio of specific heats, 7 = 5/3; Alfven and sonic Mach numbers. 
Ma = 0.25; - 3.66. 

Initial Conditions. Initially, the computational domain contains only supersonic 
and sub-Alfvenic plasma torus flow with a homogeneous spatial distribution and a 
Maxwellian velocity distribution; the pickup ions have a weak density and spherical 
spatial distribution. The magnetic and electric fields are B = Bq and E = — Uq x Bq. 
Inside Europa the electromagnetic fields are E = and B = Bq, and the bulk velocities 
of ions and electrons are also equal to zero. Here the X - axis is directed in the 
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corotation direction, the Y - axis is directed toward Jupiter, and the Z - axis is directed 
to the north, as shown in Fig. 1. In model I, cases (a) and (b) we use a homogeneous 
magnetic field for the initial and boundary conditions (see paragraph above). In model 
II, case (c) we use an extrapolation of the magnetic field profile along the E4 trajectory 
(see, Kivelson et al., 1999; 2009) onto the computation domain for the initial and 
boundary conditions. The effect of global variation on the magnetic field in the rest of 
Europa was not taken into account directly in the modeling but it was included in the 
modeling as an internal magnetic dipole (see. Schilling et al., 2007; 2008). 

At t > we begin to inject the pickup ions with a spatial distribution according 
to Eq. (8). Far upstream {x = —15 Re), the background ion flux is assumed to have a 
Maxwellian velocity distribution. 

Boundary Conditions. On the side boundaries {y — ±DY/2 and z — ±DZ/2), 
periodic boundary conditions were imposed for incoming flow particles. The pickup 
ions exit the computational domain when they intersect the side boundary surfaces 
y = DY/2 -5xAy,y= -DY/2 + 5xAy,z = DZ/2 - xAz, z = -DZ/2 + 5 x Az. 
Thus there is no influx of pickup ions at the side boundaries. 

At the side boundaries we also use a damping boundary condition for the 
electromagnetic fleld (see e.g., Lipatov and Combi, 2006; Umeda, Omura and 
Matsumoto, 2001). This procedure allows us to reduce outcoming electromagnetic 
perturbations, which may be reflected at the boundaries. 

Far downstream {x — 12 Re), we adopted a free escape condition for particles and 
the "Sommerfeld" radiation condition for the magnetic fleld (see e.g., Tikhonov and 
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Samarskii, 1963) and a free escape condition for particles with re-entry of a portion of 
the particles from the outflow plasma. 

At Europa's surface, r = Re ~ 1560 km, the particles are absorbed. In model I, 
there is no boundary condition at Europa's surface for the electromagnetic field; we also 
use our equations for the electromagnetic field, (Eqs. (2), (4) and (9) from Lipatov and 
Combi (2006)) inside Europa but using the low internal conductivity (Reynolds number. 
Re — 0.5) and a very small value for the bulk velocity that is calculated from the 
particles. In model II, we also use an inductive magnetic dipole (0, 0, —72.5) nTi?|; for 
the boundary condition at Europa's surface that simulates the effect of a nonstationary 
Jovian magnetic field at the position of Europa. In this way the jump in the electric 
field is due to the variation of the value of the conductivity and bulk velocity across 
Europa's surface. (Note that the center of Europa is at x — 0,y — 0, z — 0). 

The three-dimensional computational domain has dimensions DX = 27 Re, 
DY = 30 Re and DZ = 30 Re- We used mesh of 301 x 301 x 271 grid points, and 
5 X 10*^ and 5 x 10® particles for ions and pickup ions, respectively, for a homogeneous 
mesh computation. The particle time step Atp and the electromagnetic field time 
step AtEB satisfy the following condition: VmaxAtp < min(Aa;, Ay, A2;)/8 and 
VmaxAtEB < min(Aa;, A|/, A2;)/256. 

The global physics in Europa's environment is controlled by a set of dimensionless 
independent parameters such as Ma, A, /9e, M^/Mp, ion production and charge 
exchange rates, diffusion lengths, and the ion gyroradius e = Pci/Re- Here 
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Pci = Uo/{eB/M[c) — M\c/ujpi and the ion plasma frequency a;pi = Ai^riQe^ / Mi. Mi 
and Mp denote the ion and proton masses. For real values of the magnetic field, the 
value of the ion gyroradius is about 80 km, which is calculated from the local bulk 
velocity. The dimensionless ion gyroradius and grid spacing have the values e = 0.05 
and /^x/Re = 0.1. 

In order to study ion kinetic effects (e.g. excitation of low- frequency oscillations 
{uj « fib) by mass loading), we must satisfy the condition A < (10 — 20)c/ujph, where 
fib and cupb denote the gyrofrequency and plasma frequency for upstream ions (Winske 
et al., 1985). The above estimation of the plasma parameters shows that we have good 
resolution for the low- frequency waves (see also Lipatov et al, 2012). 

There is another problem - numerical resolution of the gyroradius on the spatial 
grid. This becomes very important near Europa's surface where the MHD model cannot 
to be used and we have to use a kinetic model to study the trajectory of heavy ions 
and their interaction with the surface of Europa. Our current model still does resolve 
this last effect and we expect to improve the model by use of a spherical system of 
coordinates in future research. 

3. Results of Europa's Environment Simulation 

3.1 Effects of plasma betas on the plasma wake structure 

In order to study the effect of plasma parameters on the structure of the plasma 

wake and the Alfven wing, we have performed modeling (model I) for two cases (a) and 
(b) with different values of the upstream background ion temperatures, pickup electron 
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temperatures, and a value of the pickup production rate near the surface of Europa. 

The following plasma parameters are chosen the same for both models: full 
magnetosphere corotation speed is Uq — 105km/s; upstream densities are po++ ~ 
10 cm-^ PS++ = 10 cm-^; magnetic field is Bq = 469 nT; B = (77.6, -140.7, -441.3) nT; 
Alfvenic Mach number Ma = 0.25; magnetosonic Mach number Mg = 3.66. The model 
of O2 atmosphere was taken from Cassidy et al. (2007), Shematovich et al. (2005) 
and Smyth and Marconi (2006). In model I, cases (a) and (b), Europa's interior is 
represented as low conducting body with Reynolds number Re = 0.5. 

Model I, case (a): upstream ion temperatures are To++ — 25 eV; T5++ = 25 eV and 
upstream electron temperature is T^^ = 20 eV. Temperatures of electrons connected with 
non-thermal and thermal pickup ions are T^^non-thermai = 20 eV; T^,thermai = 20 eV. 

Model I, case (b) (reduced density for thermal O2 by a factor 60 near surface and 
higher electron temperatures; increased upstream ion temperatures, To++ = 100 eV; 
T5++ = 100 eV): the upstream electron temperature is Tg.o — 20 eV; temperatures of 
electrons connected with non-thermal and thermal O2 pickup ions Te^non-thermai = 
200eV;Te,t;,erma; = 200eV. 

We have computed several hybrid models with different ion and electron betas, and 
different production rate for O2 pickup ions, but we discuss here only the models that 
fit the observations. 

The initial thermal velocities of O2 non-thermal and thermal ions are chosen as 
the following: Vth,non-thermai = 3.0km/s (2eV) and Vth,thermai = 0.5km/s (0.05 eV). 
The initial bulk velocity of O2 pickup ions is about Ikm/s. Eq. 8 gives the 
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following total pickup ion production rate: Qo+, thermal — 0.825 x 10^^ and 

Qo+, non-thermal — ^'^^ ^ 10^^ S ^ 

Let us consider first the global picture of the interaction of the plasma torus with 
Europa. The results of this modehng are shown in Figs. 2, 3, and 4. Figures 2 and 3 
demonstrate 2D cuts for non-thermal and thermal O2 pickup ion density profiles. One 
can observe the asymmetrical distribution of the pickup ion density (top, case (a)) and 
(bottom, case (b)) in the x-y, y-z{x — 5 Re) and z-x planes. The pickup ion motion is 
determined mainly by the electromagnetic drift. The motion along the magnetic field 
is due to the thermal velocity and the gradient of the electron pressure. A more wider 
density profile of the pickup ions was observed in the case (b). Figs. 2 and 3 (bottom). 

The figures demonstrate a strong structuring in the non-thermal and thermal O2 
ion density profiles. While case (a) produces a much higher peak in the thermal O2 ion 
density as was seen in E4 observations, case (b) produces much better agreement with 
observation for the thermal O2 ion density as shown in Figs. 2 and 3. 

The modeling also demonstrates the asymmetrical distribution of the background 

ion density in the x-y, y-z{x — 5 Re) and z-x planes. Fig. 4. The asymmetrical 
distribution of the background ions in the x-y plane may be explained by the existence 
of a strong component in the upstream magnetic field. One can also see an increase 
in the plasma density near Europa due to the formation of a magnetic barrier (not 
shown here). In case (b) this effect is stronger than in case (a). The density profiles for 
SO'^^ background ions are close to the density profiles for ions. 

The inclination of the magnetic field results in an asymmetrical boundary 
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condition for ion dynamics (penetration and reflection) in Europa's ionosphere and an 
asymmetrical Alfven wing. 

Note that the background ion flow around the effective obstacle that is produced 
by pickup ions and the ionosphere. The pickup ions flow from the "corona" across the 
magnetic field due to electromagnetic drift, whereas the motion along the magnetic field 
is determined by the thermal velocity of ions and the electron pressure. 

Figure 5 demonstrates the ID cuts {y — 0, z — 0) of the background density O"*""*" 
for case (a) (top) and case (b) (bottom). Strong jumps in the plasma density with 
No++,max — 80cm ~^ (case (a)) and No++^max — 17cm (case (b)) are observed on the 
day-side of the ionosphere, whereas a reduction in the plasma density is observed in the 
plasma wake. Note that the jump in the plasma density profile is stronger in case (a) 
than it is in case (b). Both jumps are located near the surface of Europa. 

Figures 6 shows 1-D density profiles of the background and pickup ions along the 
E4 trajectory of the Gahleo spacecraft. One can see a strong plasma void in the center 
of the plasma wake. There is also a sharp boundary with an overshoot in the density 
profiles on the side of the plasma wake in the Jupiter-direction, and a smooth boundary 
layer on the side in the anti- Jupiter direction. Fig. 6 (top). The density profile for O"*""*" 
is similar the density profiles for the S^^ upstream ions. Fig. 6 (middle and bottom) 
also shows the density profiles for the non-thermal (top) and thermal (bottom) 
pickup ions. One can see the split structure of the plasma tail. The effect of splitting 
of the plasma tail was also observed in the hybrid modeling of weak comets (see, e.g., 
Lipatov, Sauer and Baumgatel, 1997; Lipatov, 2002). The general feature of this plasma 
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density is due to the effect of the finite heavy gyroradius. The total ion density profile 
observed in E4 pass is shown in Fig. 6 (bottom). The observed value of the density 
in these peaks is lower than in modeling and it may be explained by an overestimated 
density of O2 pickup ions for case (a). In the case (b), disagreement is not as strong, an 
improvement of the atmosphere model is still required. 

The modehng gives the following total fiuxes for the O2 pickup ions (case 
(a)): 1.4 X lO^^mol/s (non-thermal) and 1.75 x lO^^mol/s (thermal); (case (b)): 
0.8 X lO^^mol/s (non-thermal) and 1.0 x lO^^mol/s (thermal) across the back boundary 
X = 12Re. 

Let us consider a global distribution of the electric and magnetic field in Europa's 
environment. Figure 7 shows B^, magnetic and Ey electric field profiles for case (a) 
(left) and case (b) (right). The y — z cuts (top and middle) are located at x/Re — 7, 
and X — y cuts (bottom) are located at |/ = 0. The figure demonstrates perturbations 
in the magnetic B^ and electric Ey field profiles, which are due to the formation of 
an Alfven wing. The increase in the magnetic field Bz indicates the formation of an 
asymmetrical magnetic barrier. Fig. 7 (bottom). 

The asymmetry of the modeling distributions in B appears to be caused by the 
finite gyroradius efi'ects of incoming and pickup ions. A weak perturbation of the 
magnetic field was observed near the ionosphere of Europa: compression of the upstream 
magnetic field and decompression in the plasma wake. 

The modeling also shows the formation of an Alfven wing in the direction of 
the main magnetic field. The formation of the Alfven wing in a sub-Alfvenic fiow 
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near Europa is similar to a formation near lo, which was first studied analytically by 
Neubauer (1980). The pickup ions play an important role in the fine structure of the 
Alfven wing due to effects of mass loading. In particular, the scale of the front of the 
Alfven wing must be determined by the gyroradius of pickup ions. Unfortunately, in our 
3D hybrid kinetic simulation we cannot yet resolved these spatial scales. 

3.2 Effects of inductive Europa's magnetic field 

In the first set of models (Sect. 3.1, model I, cases (a) and (b)), we used a 

homogeneous global magnetic field as an initial condition. These models do not produce 
agreement between the simulated and observed magnetic fields. 

In the second set of modehng we take into account the gradient of the global Jovian 
magnetic field for an initial magnetic field distribution. In the paper by Kivelson, 
Khurana, Stevenson et al. (1999); Kivelson et al. (1997); Kivelson et al. (2000), it has 
been shown that the By component of the magnetospheric magnetic field has strong 
time variations at the position of Europa. In the MHD-fiuid approximation the effects 
of such magnetic field variations are estimated in SchiUing, Neubauer and Saur (2007); 
SchiUing, Neubauer and Saur (2008). The initial plasma density and bulk velocity 
distribution in our modeling were taken from the E4 flyby data (Paterson et al., 1999). 

We created the following model II, case (c) for simulation: the density for thermal 
O2 is the same as for model I, case (b), and the pickup electron temperature is lower than 
in model I, case (b). The plasma density and bulk velocity distribution in our modeling 
were taken from the E4 fiyby data (Paterson et al., 1999): full magnetosphere corotation 
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speed Uo — 105km/s; upstream densities are po++ — 10cm~^; ps++ — 10cm~^; 
upstream ion and electron temperatures, To++ — 100 eV; T5++ — 100 eV; Tgfi — 20 eV. 
The temperatures of electrons connected with non-thermal and thermal O2 pickup ions 

are Te^non-thermal — 100 cV; Te^thermal — 100 cV. 

In our hybrid kinetic modeling (model II) we use a simple magnetic dipole model 
of the induced oceanic magnetic field from the ten-hour corotation variation of the 
background Jovian magnetic field at Europa (see paragraph "Boundary Conditions", 
Sect. 2). And, finally, we fit the results of modeling to the components of the measured 
magnetic field. 

This is not yet a fully self- consistent approach but provides a first approximation. 
Also, the ocean may not be exactly a spherically symmetric conducting shell and may 
ultimately require a higher-order multipole model for the induced fields. 

Figure 8 demonstrates the 2D cuts for non-thermal and thermal pickup ion 
densities. The figure does not show any extension of the pickup ion profile in the y 
and z directions. The plasma wake is narrower in the y and z directions compared to 
that produced by model I, cases (a) and (b). The reason for this effect is the lower 
temperature of electrons connected with pickup O2 ions than in case (b), and a lower 
pickup ion production rate near the surface of Europa than in case (a). 

Figure 9 shows the distribution of the ion density in the x-y, y-z {x — 5 Re) 
and z-x planes. The narrow plasma wake may be explained by the cooler temperature of 
the electrons connected with pickup ions, resulting in a smaller polarization electric 
field that is responsible for the expansion of Europa's ionosphere. 
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One can also see an increase in the plasma density near Europa due to the formation 
of a magnetic barrier (not shown here). The density profile for S0~^~^ background ions 
is close to the density profile for O"*""*" ions as in model I, cases (a) and (b). 

Figure 10 shows a 1-D cut of the background O'^'^ density along the x- 
axis {y = 0,z = 0). One can see jump in the background plasma density with 
^o++,max — 90cm~^ (model II, case (c)) on the day-side of the ionosphere and depletion 
in the plasma density in Europa' s plasma wake. Note that the jump in the plasma 
density profile is stronger in model II, case (c) than is observed in model I, case (a). 
The jump is located near the surface of Europa, as was observed in model I, cases (a) 
and (b). 

Figures 11 shows 1-D density profiles of the background and pickup ions along 
the E4 trajectory of the Galileo spacecraft. One can see a strong plasma void in the 
center of the plasma wake. There is also a sharp boundary with an overshoot in the 
density profiles on the left side of the plasma wake, and a smooth boundary layer on 
the right side. Fig. 11 (top). The density profile for S~^~^ is similar the density profile for 

background ions. Fig. 10 (middle) shows the density profiles for non-thermal and 
thermal O2+ pickup ions. The total ion density profile observed during the E4 pass is 
shown in Fig. 11 (bottom). Again, one can see two peaks in the total ion density profile. 
However, the observed value of the density in these peaks is lower than predicted by the 
model; this may be explained by an overestimated density of O2 pickup ions for model 
II, case (c). 

The modeling shows that the shape of Europa's global plasma environment depends 
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on a combination of the upstream plasma parameters and pickup ion and electron 
parameters. For example, reducing in the temperature of electrons connected with 
pickup ions results in a higher density of thermal pickup ions at the trajectory of 
the spacecraft (compare Fig. 6 (right) and Fig. 11). This effect is connected with 
the polarization electric field which is proportional to the gradient of the electron 
pressure. Reducing the temperature of the background upstream ions results in the 
widening of the plasma wake (compare Fig. 6 (left and right, top) and Fig. 11 (top)). 
These effects were earlier demonstrated in the 3-D hybrid simulation of Jo's plasma 
environment (Lipatov and Combi, 2006). We have found the similarities between the 
plasma environments of these objects. Indeed, lo and Europa have sufficiently thin 
exospheres and strong magnetic fields resulting in a small value of the ion gyroradius. 

Let us consider the global distribution for the electromagnetic field of model II, 
case (c). Figure 12 shows 2-D cuts for the magnetic B^, and electric Ey field profiles. 
The distributions for the B^, Ey field shown in the figure are close to the distributions 
for model I, case (b). However, there are significant differences between the B^ profiles 
for model I, case (a) and case (b), and model II, case (c). The differences between the 
Bx profiles for cases (a) and (b). Fig. 7 (top) are due to a much higher density of the 
thermal O2+ pickup ions in the plasma wake, whereas the differences between the B^ 
profiles for cases (b) and (c) are due to the nonlinear interaction of the Alfven wing with 
the inhomogeneous Jovian magnetic field in model II, case (c). 

Figure 13 shows the magnetic field components (sohd line) B^-, By , B^, and \B\ 
along the E4 trajectory of the Gahleo spacecraft. The magnetic field components of the 
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inductive magnetic dipole that simulates the effect of the nonstationarity of the Jovian 

magnetic field are shown by a dotted line ( ). The circles (o) denote observational 

data from Kivelson et al. (1997) and the initial Jovian magnetospheric field at the 
position of Europa (+++). The simulation produces a satisfactory agreement with 
the observational data for the By magnetic field component, but not for the and 
Bz magnetic field components. A multipole model for the oceanic magnetic field 
may address this issue. We will need to improve the model of the O2 atmosphere, 
the resolution of the ion trajectory, and the gradient in the atmosphere/ionosphere 
density profiles near the surface of Europa to obtain better agreement in the B^ and B^ 
magnetic field components 

4. Conclusions 

Hybrid modeling of Europa' s plasma environment for the E4 encounter with 3 ion 
species demonstrated several features: 

• The modeling shows a strong phase mixing in the plasma wake. The plasma wake 
demonstrates the formation of time-dependent structuring in the pickup ion tails 
(see, e.g., McKenzie, Sauer, Dubinin, 2001 for a weak comet case) and the splitting 
of the pickup ion tails. The splitting of the plasma wake has the same nature as 
the splitting of the weak comet's plasma wake or the splitting of Titan's plasma 
wake. Such finite gyroradius effects were also observed in 2.5 D hybrid and bi-fluid 
modeling of a weak comet (see, e.g., Lipatov, Sauer, Baumgartel, 1997; Sauer 
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et al., 1996; 1997; Lipatov, 2002) and in 3D hybrid modeling of Titan's plasma 
environment (Lipatov et al., 2011; 2012). The further investigation of these fine 
structure needs an additional modeling with much better resolution. 

• The model shows a magnetic field barrier formation at the day-side portion of 
the ionosphere. The formation of an Alfven wing in the plane of the external 
magnetic field was also observed. Note that the Alfven wing was earlier observed 
in a hybrid simulation of the plasma environment of lo and Europa by Lipatov 
and Combi (2006) and by Lipatov et al. (2010). An MHD simulation of the 
plasma environment of lo and Europa also produces the formation of an Alfven 
wing (Saur et al., 1999; 1998; Liu et al, 2000; Schilling et al., 2008). 

• The ion and electron temperatures play an important role in plasma structure 
formation, and in creating the ion fluxes inside the ionosphere. These effects were 
observed earlier in a 3-D hybrid simulation of lo's plasma environment (Lipatov 
and Combi, 2006). The hybrid model produces the correct pickup ion flux along 
the magnetic field, in contrast to the MHD models which operate with pickup ions 
with a Maxwellian velocity distribution. In the current paper we have presented 
only three runs with different combinations of the upstream ion and pickup 
electron temperatures. 

• The model's total ion density in the plasma wake does not satisfactory match the 
observed density. 
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• The constant induced dipole moment (model II, case (c)) improves a fit of the 
magnetic field By component to the E4 trajectory. However, a fit of the magnetic 
field component is still not satisfactory due to the imperfect model of the 
atmosphere/ionosphere and unsatisfactory numerical resolution of the gyroradii 
on the grid cell. 

• Use of an inhomogeneous background magnetic field provides a good agreement 
between the observed and simulated magnetic fields. However, we still need to 
improve the resolution of the gradient in the atmosphere density, the gyroradius of 
pickup ions, and the fields in the internal non-conduction ice shell and conduction 
ocean layers of Europa. 

In our future computational models, we plan to include a nonstationary boundary 
condition for the magnetic field in order to take into account the spatially inhomogeneous 
and nonstationary background Jovian magnetic field. This model will also be appropriate 
for a potentially nonspherical ocean shell. We also plan the use of a varying atmospheric 
density, a varying electron temperature (that plays key- role in the pickup ion dynamics) , 
and sputtering processes (Johnson, 1990; Johnson et al., 1998) at the surface of Europa. 
We also plan to use a composite grid structure using the "cubed sphere" technique (see, 
e.g. Koldoba et al, 2002) to improve the resolution of the a small scales near the surface 
of Europa and to increase the size of the computational domain. 

The composite grid structure will allow us to estimate the inductive magnetic field 
from the ocean as a part of the total current closure that also includes the external 
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plasma currents. This technique will allow us to study wave-particle interaction effects 
in the far plasma wake, such as ion cyclotron waves that have been observed in the 
Galileo fiyby mission (see e.g. Volwerk, Kivelson and Khurana, 2001; Kivelson, Khurana 
and Volwerk, 2009). These models must include the induced magnetic field from a 
putative subsurface ocean, and will also include particle trajectory tracing for test 
particles, e.g. electrons and high-energy ions. 

Note that the larger computational domain allows us to use the upstream 
parameters for the plasma and electromagnetic field instead of the use of the "damping" 
boundary condition. However, in the outer region of the computational domain (large 
cell size) we have to use a drift- kinetic approach (see e.g. Lipatov et al., 2005) for ion 
dynamics since we cannot approximate the ion trajectory there. We can also use a 
complex particle kinetic technique (see e.g. Lipatov, 2012) which provides a flexible 
fluid/kinetic description and may significantly save computational resources. 
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Figure 1. Europa's environment and system of coordinates. 
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Figure 2. 2-D cuts of non-thermal pickup ion density profile. Model I, case (a) (top) and 
case (b) (bottom), x — y cuts (left column) are located at z = 0, y — z cuts are located at 
x/Re = 7, and x — z cuts (right column) are located at y = 0. 
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Figure 3. 2-D cuts of the thermal pickup ion O2 density profile. Model I, case (a) (top) and 
case (b) (bottom), x — y cuts (left column) are located at z = 0, y — z cuts are located at 
x/Re = 7, and x — z cuts (right column) are located at y = 0. 
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Figure 4. 2-D cuts of the background O^^ ion density profiles. Model I, case (a) (top) and 
case (b) (bottom), x — y cuts (left column) are located at z = 0, y — z cuts are located at 
x/Re = 7, and x — z cuts (right column) are located at y = 0. 
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Figure 5. ID cuts of the background O"*""^ ion density profile. The cuts are located at y = 0, 
z = 0. Model I, case (a) (top) and case (b) (bottom). 
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Figure 6. ID cuts of the background and pickup non-thermal/thermal (Oj) ion 

densities from simulation. Y{Re) denotes a projection of satellite trajectory onto the y axis. 
Model I, case (a) (left) and case (b) (right). Bottom - E4 observation of the total ion density 



(Paterson et al. 1999). 
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Figure 7. 2-D cuts of the B^, magnetic and Ey electric field profiles. Model I, case (a) 
(left) and case (b) (right), y — z cuts (top and middle) are located at x/Re = 7, and y — x cuts 



(bottom) are located at z/Re = 0. 
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Figure 8. 2-D cuts of the pickup ion density profile. Non-thermal (top), thermal 
(bottom). Model II, case (c). x — y cuts (left column) are located at z = 0, y — z cuts are 
located at x/Re = 7, and x — z cuts (right column) are located at y = 0. 
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Figure 9. 2-D cuts of the background O"*"*" ion density profiles. Model II, case (c). x — y cuts 
(left column) are located at z = 0, y — z cuts are located at x/Re = 7, and x — z cuts (right 
column) are located at y = 0. 
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Figure 10. ID cuts of the background O"*"*" ion density profile. The cut is located at y = 0, 



z = 0. Model II, case (c). 
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Figure 11. Background {S^^, O^^), non-thermal (Oj) and thermal (Oj) pickup ion densities 
from simulation. Y{Re) denotes a projection of the spacecraft position onto the y axis. Model 
II, case (c). Bottom - E4 observation of the total ion density (Paterson et al. 1999). 
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Figure 12. Model II, case (c). 2-D cuts of the magnetic (top), electric Ey (top) and 
magnetic Bz (bottom) field profiles, y — z cuts (top) are located at x/Re = 7, and y — x cuts 



(bottom) are located at z/Re = 0. 



49 



150 
^ 100 

" 0^ ^ 

-50 

-2 2 

100; ^ 

^ 0^ 

^ -100^ 

" -200 h 
-300 i 



200- 

^ 0- 
H 

-200- 

N 

" -400 
-600 



H 300 r 
ffl 200 h 
100^ 

oL 







2 

X(Re) 



Figure 13. Magnetic field component profiles along the E4 trajectory after fitting with induc- 
tive dipole magnetic field. Solid line - modeling, (-) denotes dipole field, and (-I-) is the Jovian 
magnetic field at the position of Europa. o - Galileo's E4 flyby measurements (Kivelson et al. 
1997). X{Re) denotes a projection of the spacecraft position onto the x-axis. 
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